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Abstract 

A numerical method to solve linear integro-difFerential equations is presented. 
This method has been used to solve the QCD Altarelli-Parisi evolution equations 
within the HI Collaboration at DESY-Hamburg. Mathematical aspects and numer- 
ical approximations are described. The precision of the method is discussed. 



This article is an extended version of an unpublished note 

In a recent publication j^J, P. Ratcliffe proposed a numerical method similar to the 
one described in ref. J^] . In addition, he pointed out the problem of non commutativity 
of the Next-to-Leading-Logarithmic-Approximation (NLLA) Altarelli-Parisi kernels and 
the fact that we did not account for in ref. f^] . Our present extension of jl^/ therefore 
concerns the account for non commutativity. 

After having given the proper modification of the numerical method (section 2.1) we 
explicitly show that non commutativity effects can safely be neglected provided the 
evolution is performed, as usual, from points to points on a grid (section 2.2). The rest 
of the paper is untouched, in particular the references are not updated. 



Introduction 

Inclusive Deep Inelastic lepton-hadron Scattering (DIS) cross section measurements offer a 
powerful test of perturbative Quantum Chromo Dynamics (pQCD) The DIS process 
i{k)'^h{P) — »• i^{k')X (here we shall only consider the case of charged leptons in order 
to simplify the discussion) kinematic is described by two Lorentz invariants. One usually 
chooses the transferred momentum squared = —q^ = {k' — k^ and the Bjorken variable 
X = Q'^/{2P.q). In terms of these two kinematic variables, the internal dynamics of the 
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struck hadron enters the cross section via three structure functions: Fi{x,Q^), F2{x,Q^) 
and F^^XjQ"^). Because only F2{x,Q'^) contributes significantly to the cross section for 

< M'^^c^, we shall concentrate on this structure function in the present paper. 

Within the framework of pQCD, F2{x,Q'^) is given by the convolution in x of the 
well known Wilson coefficients and of the parton densities inside the hadron (we shall 
work in the MS factorisation and renormalisation scheme). The densities of partons, 
consisting of quarks and gluons, are computed from the solution of the Altarelli-Parisi 
(AP) equations [|l|. According to 0|, it is useful to define the gluon g, a non-singlet q^s 
and a singlet S quark combination densities. They are the solution of the set of AP 
integro-differential equations: 

dqNS dw 



dt 



I ^P^,{w,t)qMs{^,t) (1) 

J X ^ 

(2) 



dt\g{x,t)) w \Pg,{w,t) Pgg{w,t) J \g{x/w,t) 

where all kernels P are expanded perturbatively: 

P = a,(t)pW(^) + a^(t)p[2](^) + 0(^3) (3) 

Here t = log((5^/A^) and A ^ 200MeV is the pQCD scale parameter. Expressions of 
the leading and next leading order splitting functions, P'^^ and P'^^, can be found in P|. 
The expression of the strong coupling constant from ref. [^] will be used. However, 
eq. (|lD,(@) hold only for ^ A^ where the perturbative series (^) is convergent. Some 
non-perturbative input functions are then required to solve the system. In practice, these 
functions depend on unknown parameters (except in the case of [Q) which are determined 
from a fit to experimental structure function measurements Thus, eq. (|I|),(0) must 
be solved many times during the usual minimisation procedure and a fast and precise 
numerical method is then required. 

Two [jdifferent numerical methods have been used so far to solve eq. (P,(|D (see 
ref. [§] for a review). The first one uses the fact that the Mellin transform of the 
system (|1|J^) leads to a simple set of first order differential equations in the x complex 
conjugate moment space. The solution is straightforward, after some further necessary 
simplification of the evolution kernels, but the price to pay is to perform numerically an 
inverse Mellin transform. The method proposed in the present article belongs to a second 
kind of approach in which the system is solved in the x space: let us comment with more 
details the main features of three existent approaches of this kind. 

• It is first natural to use methods based on the Taylor expansion of the parton 
densities in log(Q^) [|, |n| (called 'brutal force' method in ref [§, |n|). But it is 
well known [0 that such techniques lead to rounding errors when the evolution 
covers many order of magnitude. And most of all, reaching a good accuracy is 



prohibited by a necessary high CPU (see ref. |TT|] for a complete numerical study). 



^ Analysis of F2{x,Q^) moments will not he considered here. It is well known since a long time 
that it relies too much on the behaviour of the structure function in the elastic a; w 1 and 'small-x ' a; ^ 
regions. 



2 



It was early observed [|T3] that shifted Jacobi polynomials offer many useful 
properties to solve AP equations: orthogonahty for x G [0, 1], a weight function 
x^{l — x)" describing the asymptotic behaviour of parton densities (a > 0, /3 > — 1 
are real arbitrary numbers) To understand advantages and disadvantages of that 
method, let us write the expansion of the non singlet densities in terms of the 
shifted Jacobi polynomials 0^'^(a;): 

oo 

XQNsix, t) = - xr J2 W ^ ^n^ix) ■ 

n=0 

where 9"'^(x) = J2]j=o(^j,n{ci, P)^-' ■ Using orthogonality relation one immediately 
relates the coefficients a„(t) to the non singlet moments. Although the solution 
can be explained in a very compact form, the sum in the previous equation must 
be truncated |TH]: when taking for example the smallest experimental x = 10~^ 



accessible value |]T^, one sees that increasing n leads to numerical rounding errors. 
Note that this method was originally used for analysing data at x > 10^^ and, as 
pointed out in 0, 'a careful study of numerical stability is required to extend this 
method at lower x values'. Recently, it was claimed [1^] that this method could be 



applied down to x ~ 10 but no numerical stability studies have been yet provided. 



In ref. |T^, the authors made a series expansion of both parton densities and kernels. 



Changing the integration variable x to log(l/x) they show that the best polynomial 
basis to perform this expansion are the Laguerre orthogonal Polynomials |T^. In 
particular, the convolution product property of these polynomials leads to an expli- 
cate mathematical form of the solution. Using the generalised Laguerre polynomials, 
they also include the asymptotic behaviour of the parton densities as depicted in 
the previous item. The restrictions of this method are twofold: the kernels and the 
parton densities are approximated leading to series in the solution which must be 
truncated. Furthermore, some recent studies of the numerical precision |TI], M have 



shown that this method 'may not be satisfactory at small-x and large-x' |11 



With the forthcoming high precision F2 measurements by the HERA experiments p6| , 
together with the very accurate fixed target experiment data[^ |2l| , precise solutions 



of AP equations are required in a kinematic plane covering five orders of magnitude in 
X and Q^. In order to perform QCD analysis in this large kinematic domain, we present 
a new numerical method to solve the system (|I],0). Using the finite element method, 
together with a scaling property of the convolution integrals appearing in eq. (0),(|D, we 
show that it is possible to obtain an explicit formula for the solution of AP equations 
without any infinite series. We also show that our formulation leads fortunately to a 
reduction of the CPU time by one order of magnitude with respect to the 'brutal' method 



described in [11 



The rest of this article in organised as follows. The method is described in sect. [l|: 
notations and formalism of are used. Application to the AP evolution equations is 
described in sect. ^. The numerical precision, together with the performances are studied 
in sect. 131. 
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1 Description of the method 



Let us write generically AP integro-differential equations (|T]),(0) in the form: 

dF{x,t) 



dt 



(4) 



where F stands for quark or gluon densities in the proton and therefore belongs to a class 
of smooth functions. In addition, one has x G [0, 1], t G [0, oo] and = for all t. 

We start to solve eq. @) by defining an ordered sequence of points >= {xi}"^ and 
the corresponding sequence \F >= {Fi}^ with Fi = F{xi,t) . It is then possible to define 
a piecewise polynomial function ^r{x) of order r which interpolates |F > in x at a given 
value of t: dri^i) = Fi. 

^'^ is described by a set of polynomial functions |^ >= of order r: 

dr{x) = ^[(x) if Xi <X < Xi+i. 

together with \x > define a linear space Wr^x- If one supplies a sequence of continuity 
conditions \u >= {i^i}2~^^ the defined space 'iPr,x,u is still linear and is a subspace of 
Pr,a;- Because of the linearity of ^r,x, one can define a basis of continuous functions 
10 >= {0j(x)}" such that (pi{xj) = 6ij stands for the kroeneker symbol). In this way, 
one can write 

n 

^r{x) = J2F,Ux), (5) 

i=l 

with |0i >G fr,x,u- Equation is valid if and only if 

dim fr,x,u = r X n — Vj = n , (6) 

i=2 

and is thus only satisfied for a restrictive class of sequences \v > . Functions '^r which are 
constructed in such a way are the original spline functions . 
We turn now to the solution of (|D by making two assumptions: 

/■I dw 

'^r{x) = F{x) and / — K{x/w)(j)i{w) G '^r,xM- 

Jx W 

Depending on the choice of |x > and r, theorems concerning the precision of these ap- 
proximations are given in p3[. This topic will be covered in sect. |^ in the case of AP 
equations. 

Using (H), eq. takes a matrix form 

d ^ 

■F>=Y,<{t)Mm\F>, (7) 



dt 



m=l 



'We use the following convention: {a^i}" = {xi,...,x„ 



4 



with 



< t\M^\j >= I ' —pl-^^{x,/w)<j),{w) 



„. w 



The solution of (J^) can therefore be written formally 

|F>=exp[^ t dt'aT{t')M,^\Fo>, (9) 

where \F >= \Fq > at t = to {to > can be chosen arbitrary). 

Next we choose the sequence \x > and the basis functions |0 > according to the 
condition (^. We also focus on the numerical computation aspects, namely the CPU 
time and the precision. 

We first remark that if 



= Xxi, with A > 1 and 




0i+i(Ax), 

for a; > Xj+i, 



matrices Mm are triangular and < i + l\Mm\j + 1 >=< i\Mm\j >■ This very interesting 
scaling property relies on the fact that the integrations are convolution products running 
from X to 1; it remains true even in the presence of a '+' (QCD) prescription appearing 
in the AP kernels [0. 

As we are free to choose the order r of the piecewise polynomials, we take the simplest 
case, i.e linear interpolation (r = 2). It is known to be nearly optimal for large n values. 
The corresponding basis is unique and consists on the Lagrange (or sometime called 'hat') 
functions: 

{{X - Xi_i)/{Xi - Xi^i) , Xi_i <X<Xi 
{Xi+i - x)/{Xi+i - Xi) , Xi<X < Xi+i (10) 
, otherwise. 

Note that integrals of eq. (^ are simple to compute using definition (|TD|); cancellations 
of singularities contained in the AP kernels can even be computed analytically. 

If one wishes to design a computer program it is necessary to define completely the 
points where the function F has to be calculated. This is done by defining a sequence in 
the variable t\ \t >= The solution may then be propagated from the beginning of 

the grid to its end by using a recursive expression in analogy with eq. (1): 

\Fk >= expAk\Fk-i > with Ak = J2 dt'a'^{t')Mm , A; = 1, . . . ,nt . (11) 

m=l 

In the case of a fitting procedure, if the kernel K{w,t) is not changed [i.e as fixed at a 
certain scale t) it is sufficient to compute the matrices exp{Ak) only once. Furthermore, 
the type of integrals appearing in eq. S are computed at the initialisation step because 
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they depend on the choice of |x > only. They can be computed with a high numerical 
precision. If only |Fo > is changed one must only redo the matrix multiplications implied 
by eq. (0): this is a consequence of the linearity of (^). 

To efficiently exponentiate matrices Ak, it is powerful to introduce the 'band' matrices 
B defined by: 

< i\Bi\j >= 6i+ij , />0 . 
It is easy to show that B matrices fulfil the multiplication rule 

BiBj = BjBi = Bi^j. 

These matrices form a basis for the kernels matrices and we can write 

n-1 2 tk 

Ak = J2 ^(^iBi with ^ai = Y^ / dt'a'^it') < /|M^|1 > . 

1=0 m=l 

The product of two matrices such as Ak will be a matrix of the same type, z.e a linear 
combination of the i?'s. The use of band matrices leads to a number of operations of 
order whereas ordinary matrices would give a CPU time increasing as n^. 
Let us now isolate the diagonal term in A^ and rewrite: 

n-1 

Ak = J2^aiBi = ''aoBo + Al 

i=0 

Commutativity allows us to split the exponential in two parts: exp(Afc) = expC^ao) exp(y4'^). 
In this expression the second exponential can be expanded as a series which will consists 
in a sum of products of B matrices. As Bq is not an element of those products and be- 
cause of the multiplication rule the index of the 5's will increase along the sum to reach 
its maximum n — 1. So there will be a limited number of terms in the series and one can 
write 

n-1 (A/ Y n-1 

exp(A,) = expC^ao) E = E '^^^r (12) 
1=0 ■ j=0 

If we turn back to eq. (0), we can finally write the solution in the following form: 

n—i 

F{xi, tk+i) = J2 tk). (13) 

j=0 

This equation shows explicitly that F{xi, Q"^) depends only on the values of the function 
F for X > Xi. 
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2 Application to the Altarelli-Parisi equations 

It is clear that the non-singlet equation ( pTSD is a solution of eq. (|l]). In the case of singlet 
and gluon coupled differential equations, we may still retain formally eq. (|^) but A is now 
made of four sub matrices. Let us rewrite eq. (H) with slightly different notations: 

>= e-^'\Tk > with A = f^-^^ ^"^"j , (14) 

where the ^matrix can longer be expended linearly on the basis of band matrices, unlike 
Aqq, Aqg, Agg and Agg. We can also consider that the vectorial space involved is a direct 
product of a n-space and of a 2-space and write (see Appendix A) : 

A = AqS + ^ = AqS + A^X + Ay-y + A,Z , 

where 

The following rules hold: 

£ = £-i = X'^ = -y'^ = , (15) 
Xy = Z + permutations , (16) 

(17) 

S commutes with X, y, Z and those anti commute between themselves. Therefore, one 
can split into pieces the exponential 

= e^°e^ . 

The second exponential in the r.h.s of the above equation is computed by series expansion. 
Then, one needs to calculate powers of A . Making use of eq (|T3p one gets (see Appendix 
A): 

{1y = [A^X + Ayy + A^Zf = [Al - Al + A^^£ = A\E. 

This shows that we can split the series in its even and odd parts and sum them indepen- 
dently. The result is the following: 

= f ^4^^ + '-^^A-^-A 1 . (18) 

y4g, its square root A^. and the inverse are easily computed fusing the i?'s multiplication 
rules. 



^Depending on the sign of 's diagonal, will be real or pure imaginary; it turns out for AP 
equations that we are in the first case. 
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2.1 Non commutativity correction 



Eq. is a solution of the AP equations if the following relation holds: 

d A d . A ., A 
dt dt 

Looking at the series expansion of the exponential and its derivative one can see that it 
is the case only at the leading order in where -^A and A commute. 

In the general case, i.e. including higher orders is a^, the same form for the solution 
may be kept by defining an operator B such that 

^"•^ > - > and \T >= e^l^ >U ^e^ = A' x 



dt ' ' ' j dt 

Simple but tedious non commutative algebra (see Appendix A) gives: 



where H is given by eq. As those correction terms are always small (see below) 

integration can be made by a simple trapezoidal formula onto a grid and the system 
can be solved iteratively using A! , A as initial value for B' , B. 



2.2 Effect of the non commutativity correction 

To check the validity of this approach the singlet density have been evolved in the NLLA 
according to four different assumptions: 

1. the solution is propagated from one point of a grid to the next one as usually 
done and the non commutativity correction is not applied. 

2. the solution is propagated from one point of the grid to the next one and the 
non commutativity correction is applied. This method should be the best one and 
is used has the basis of the comparison. 

3. the solution is propagated from the first point of the grid to the actual one and 
the non commutativity correction is not applied. No grid is considered. 

4. the solution is propagated from the first point of the grid to the actual one and 
the non commutativity correction is applied. No grid is considered. 

To compare the different results we perform the ratio of the singlet densities ^2(2;, Q^)/Ei(x, Q^), 
^3(0;, (5^)/Si(x, Q^) and T^i{x,Q'^) /Tii[x,Q'^) where the subscript refers to the above as- 
sumptions. As for the gluon density, the effect is found to be an order of magnitude 
smaller than for the singlet case. 

Figure |l] shows the singlet ratios as function of x at = 30000 GeV^ and as function 
of at X = 0.0001 (S2(x,g2)/Si(x,Q2) is scaled by a factor 1000 in this figure). The 
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Figure 1: 


Non commutative corrections, 


see text. 



starting scale of the evolution is 1 GeV^. The parametrisation of the gluon and singlet 



densities at this scale is taken from a recent global fit |27 . 

As expected method 1 and 2 are very close, the bias being less than 10~^ with an 
irregular shape showing that it comes entirely from rounding errors. For method 3 the 
bias is noticible at small x and large where it goes up to 6.5 %. When the non 
commutativity correction is applied (method 4 ) the bias is under control and less than 1 



This study shows that, if one propagates the Altarelli-Parisi solutions from 
point to point on a grid, the non commutativity bias is completely negligible 
and that it is not even necessary to apply the correction. 



^ This remaining bias is independent of the x and grid sizes and seems due to computer accuracy 
(the number of convolution products needed being realitvely high) 
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3 Precision and performances 



The precision of the method proposed in this paper can be easily estimated by changing 
the number n of breakpoints in \x >. The convergence of the evolved parton densities 
leads to an estimation of the precision. Input functions and A (A'^ = 225MeV) value of 



24 are used: 



xg{x) = 1.86x-°-22(l 

= 1.15x"°-^^(l - + 3.12x) 

XQNsix) = 1.14a;°-^^(l - x)^-^^{l + 8.68a;) 

at Qo = 4 GeV^. Here, T, = u + u + d + d + s + s and qns = u + u — T,/3 is the chosen 
non-singlet quark density (for more details concerning the fit we refer to [p4[ ). 

We define the two bounded sequences, \x > such that Xi G [xmin = 10~^,1],« = 
1, . . . ,n , and \t > such that tj G [0,10],j = 1,... ,nt. We start from n = 100 and 
increase up to n = 900 by steps of 100. rif is set to 100 for all n and this does not 
introduce numerical uncertainties since the expression of eq. ( ]T^ ) is exact. Evolutions of 
eq. (|l|),(0) are performed up to = lO^GeV^. Fig. (||a) shows the ratio Rq^^{x^n) 
defined by: 

Jx, n) = ^^^yi"'"^ , for n = 800, 700, 600, 500, 400, 300, 200, 100 (19) 

qNS[X)\n=900,nt 

for = lO'^GeV'^. Note that the calculations done with large values of n are stable 
because of a choice of a first order spline interpolations (see section |1]) . With higher order 
splines, the choice of the breakpoint set \x > is very important and solutions may be 
unstable |]23|. Figs. (0b, c) show R-s{x,n) and Rg{x,n) respectively. These plots illustrate 
the convergence of the method. The wiggles appearing in these plots are generated by 
the Lagrange functions which are also used to interpolate between the breakpoints of 
the sequence \x >. From these plots, one sees that for n = 300, the precision is of the 
order of 0.5% in the 'low-x' region. This is accurate enough comparing to measurement 



uncertainties |]T6[,|2y]. However, at 'high x\ huge differences appear. Although these 
instabilities of the calculations do not modify the results in the 'low-x' region, one cannot 
use this method to compute the structure functions at 'low-x' and 'high x' simultaneously 
with a good precision. A straightforward modification of the method is to define a second 
net in x. This new net starts from = xH^i- Then n' breakpoints are set equidistantly 
in log(x) between two adjacent points of the first net up to x = 1. For the calculations, 
the second net is first used to solve the AP equations. The extension from down 
to Xmin is then made by solving AP equations using the first net. For example, taking 
Xmin = 10^^ and n' = 5, the precision for x > = 10^^ will be equivalent to the use 
of one net made of n x n' breakpoints. Note that the CPU time is only multiplied by two 
instead of {n')"^. To show the resulting precision the fit described above is performed with 
n = 700 and n' = 1, 4, 8 {x'^^^ = 10'^ 5.6 x 10"^ 2.4 x 10"^). Then the ratio 

R' n' n) = 'l^six)\n',n,n, 
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is computed. For n = 700 and n' = 8, one can reasonably assume that the precision is 
nearly optimum at x > 2.4 x 10~^. We shall take the parton densities, computed with 
this conditions, as the reference 'high-x' sets. Fig.(^) show Rq^^g, R't. ^"^^ R'g = '''OO- 
the gain in precision at 'high- a; ' is very significant when going from n' = 1 io n' = 4. 
On Fig.(0a,b,c), the ratio R' for n = 250 and n' = 1,4,8 is also shown. As pointed out 
above, n = 250 leads to structure function computations precise enough at 'small-x'. One 
can then observe that setting n' = 8 leads to great improvement of these computation at 
'high-x'. We conclude that for n = 250, n' > 4, it is possible to perform a very precise 
calculation of the structure functions within five orders of magnitude in x and we point 
out that any global pQCD structure function analysis should pay much attention to the 
numerical precision (as it was already mentioned in [§). 

One can design a very fast (CPU) and efficient procedure to determine the input 
functions from a minimisation: first the minimum is approached rapidly by setting 
n = 100, then the 'true' minimum is reached by increasing n in steps of 100 until n = 250. 
If 'high x' measurements enter the fit, a second net must be considered at the end of this 
procedure. 

According to the AP equations, the momentum sum-rule 

S{t) = J x^{x,t) + g{x,t)^dx, 

should stay constant during the evolution in t. This is a sensitive test of the computational 
precision. With the parametrisations chosen for the input functions [|^, S{Ql) is com- 
putable analytically and is equal to a sum of Beta functions. S{Ql) is then set to 1 and 
S'((5^) is computed numerically with a linear extrapolation in the interval x G [0, 10~^]. 
This extrapolation introduces a small bias in the calculation of S'(Q^). Fig. ^ shows 
S{Q^) as a function of for the nine sequences of breakpoints of eq. (pJ]) (here one x 
net is considered). In any case, the deviation from 1 for an evolution in t which covers 
four orders of magnitude is of the order of 2/10000. 



Conclusion 

We have presented a numerical method to solve the Altarelli-Parisi equations in x space. 
Unlike conventional methods based on Taylor expansion, the evolution is computed 
'exactly'. All convolution products involving Altarelli-Parisi Kernels are computed only 
once with high accuracy. Using the Lagrange functions to construct a basis and taking 
advantage of scaling properties of AP equations, the number of operations increases with 
instead of like for other simple basis functions. In some senses, our method is 



equivalent to the polynomial approaches and specially the one of ref. |]T8[ . But, the x space 
discretization allows the use of first order polynomials to avoid precision computational 
problems and the freedom in the choice of the breakpoint sequence leads to a formal 
solution involving only finite series. The only numerical instability comes from the number 
of breakpoints in the \x > sequence. The resulting numerical precision is found to be very 
good when a large number of breakpoints are considered or when two \x > sequences are 
used. 
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A Appendix 
A.l Notations 

We shall adopt ad hoc, non standard, notations. All the components Ai and Bi with 
i — 0,x,y, z can be expended linearly on the basis of band matrices. 

• A = AoS + ^ = AoS + A^X + Ay-y + A,Z] 

• Explicitly: 



/ . 


. 


^0,x ■ 


Afi^x^ 


. 


. 


. 


■ Aq^x 




A 


. 


. 


\ . 


■ ^0,x 


. 


■ ) 



where Aj^-c are the components of the matrix A^ (there is only n independent values) 
and n + 1 is the number of discrete points in the x space. 

• A = AxX + Ayy + AzZ is an element of the three dimensional vectorial space 
generated by X, y,Z. 

. [A,B] = [1,^]. 

• {~^^)o — {Ar^Br^ — AyBy + AzBz)£ Is the £ component of the product of ~^ and 

• is the vector component of the product of ~^ and ^ . 

• A^ ^J(fl\ = a/T^ = AeS with Al^Al-Al + Al 

(A^ —AxAy AxAz\ 
AxAy —Ay AyAz is a matrix in the vectorial space. 
AxAz —AxAy A?^ j 

A. 2 Theorems 

It is easy to demonstrate the following. 

2. ^-W^ 
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4. {A B B)o = 

5. 0^)0^ = (1 

6. {1 A^)^ = A^) 

7. (l + ll A ~^)-^ = 1-1^ A '^(1 + lA^)'^ 



A. 3 Solution of f e'^ = A' x e 



Using 

^ =A' xe 



dt 

and the commutation relations of the Ai band matrices on gets 

(e^y = (e^-^°)' = (^'-5^)e^. 

We will now assume that B'q = Aq and show later its justification. As a consequence this 
equation becomes: 

(e^)' = (21) 

In order to use eq. 0, let us define: 



G 



B _ -B 

F= B-\ 

2 

H = GP-^ - 1 , (22) 
so that G' = FBB' and F' = FHB'B-\ Equation ^ becomes: 

FBB' + FHB'B-^lf + F^' = ~1'{G + F^) 

or equivalently 

BB' + HB'B-^1^ + 1^' = ^'{1 + H + ^). (23) 
The component reads 

BB'=0^')o, (24) 

and the vector component 



HBB'B-^lf + 1^' = {1 + H)'l' + '1'^ . (25) 
13 



Multiplying the two sides of this equation by B and taking the component one gets: 



(1 + H)BB' = (1 + H)(l'^\ + , (26) 

using 

BB' = \{B^y = \{n^)', = {n^\. 

The last term of eq. ^ being nul (theorem 4) we find the result of eq. ^ (i.e. the 
component). Hence, among the four equations only three are independent: having fixed 
B^ = A'qWQ obtain a system of three equations with three unknowns (the components of 
B '). Therefore, in virtue of the unicity of the solution of a first order differential equation, 
B'q = A'q is justified a posteriori. 



Equation |25| may be rewritten 

(^1 + HB-^1^ A ^' = (1 + H)^' + '1'^ 
and solved with theorem 7 

^' = - HB-^^ A ^(1 + H)-^^ (^{1 + H)'l' + 1"'^^ 
A further simplification, using theorem 4, leads to the final expression: 



(27) 



B' - A' = H A' + A' B - HB^'^iB A')oB . (28) 
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Figure 2: a) _Rg^^(x,n); b) R^{x,n); c) Rg{x,n). See eq. ( |7^ /or definitions. The less 
accurate curve (far from the line R = 1) corresponds to n = 100. The next one to n = 200 
etc... 
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Figure 3: a) Rg^^{x,n,n') ; b) R'j.{x,n,n') ; c) R'g{x,n,n'). See eq. for definitions. 

The full lines correspond to n = 700 and the dashed lines correspond to n = 250. The less 
accurate curve (far from the line R = 1) corresponds to n' = 1. The next one to n' = 4. 
The last one (defined only in the case n = 250) to n' = 8. 
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Figure 4: Momentum sum-rule as function ofQ^. The curves described in the teoct super- 
pose almost exactly. 
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